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^+ ■ The relationship between spatially heterogeneous dynamics (SHD) and jamming is studied in 

f^ ' a glass-forming binary Lennard- Jones system via molecular dynamics simulations. It has been 

f^ ' suggested by O'Hern et al. [J| that the probability distribution of interparticle forces P{F) develops a 

fS) , peak at the glass transition temperature Tg, and that the large force inhomogeneities, responsible for 

structural arrest in granular materials, are related to dynamical heterogeneities in supercooled liquids 

that form glasses. It has been further suggested that "force chains" present in granular materials 

may exist in supercooled liquids, and may provide an order parameter for the glass transition. Our 

goal is to investigate the extent to which the forces experienced by particles in a glass-forming liquid 

are related to SHD, and compare these forces to those observed in granular materials and other 

glass-forming systems. Our results are summarized as follows. We find no peak in P{F) at any 

temperature in our system, even below Tg. We also find that particles that have been localized for a 

P^ ' long time are less likely to experience high relative force and that mobile particles experience higher 

P5 ^ relative forces at shorter time scales, indicating a correlation between pairwise forces and particle 

mobility. We construct force chains based on the magnitude of pairwise forces. We find that force 

chains constructed in this manner are composed of both localized and mobile particles, therefore 

there is no one-to-one correspondence between force chains as defined here and locally mobile or 

immobile regions of the liquid. We also find that force chains do not play the same role as force 

Cy ■ chains in granular materials, but may indicate a difference in the evolution of the local environment 

of particles with difi^erent mobility. We also discuss a possible relationship between force chains 

found here and the development of string-like motion found in other glass- forming liquids 0, 3 • 

?H ' PACS numbers: PACS numbers: 64.70.pf, 61.20.Ic 
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Attempts are underway in the statistical mechanics community to unify concepts regarding supercooled liquids and 
granular materials, two very different classes of systems that display similar behavior in many respects. Systems near 
^1^ ] their glass transition, colloidal suspensions at large pressures or densities, foams under shear, granular materials under 
\^ • "tapping" or shearing, and other systems like bubbles and droplets "jam" under certain circumstances. Jamming is 
^^ ■ also seen in polymer crazes ,4]. Jamming is a process in which systems appear "stuck" in phase space because their 
\l , particles come in close contact with each other resulting in structural arrest. This process sketches what happens 
^^ ■ with supercooled liquids at the glass transition temperature Tg, where jamming is controlled by the temperature or 
density. In colloidal suspensions under high pressure or at high density, particles also get "stuck" , and the suspension 
acts as an amorphous solid 5|. Foams or emulsions flow under high shear rate, but at low shear they stop flowing and 
appear to be solids as wellQ- Extensive reviews of granular materials can be found in Refs. B, B, B ■ A mong many 
experiments that deal with the behavior of granular materials under external perturbation, Refs. |l0lllli ri2| indicate 
(~| ■ that the response of granular materials upon isolated tapping or continuous vibration is similar to the response of 
O ■ supercooled liquids near Tg. In these experiments, the density fluctuations of a granular material subject to shaking 
O ■ were investigated by reducing the shaking frequency until the granular material reached a "jammed" configuration. 
^ ' What are the similarities and differences between glass-forming liquids and granular materials? First we note some 

of the properties of granular materials. Unlike glass-forming liquids, granular materials consist of a large number of 
particles that are individually solid (grains). The grain-grain interactions are classical because the size of the grains 
^ is much larger than the de Broglie wavelength. The grains exert forces only when they are in contact, and may be 
surrounded by a fluid (typically air) or a vacuum. The collisions between grains are in general inelastic. The main 
difference between liquids and granular materials, aside from the obvious difference of length and time scale, is the 
concomitant difference in energy scales. In granular materials, thermal energies ki,T are insignificant compared to 
the energy it takes to move a single particle. Thermal energy in the liquid allows it to explore different states, while 
the low relative magnitude of thermal energy in a granular material does not allow it to sample other configurations 
unless energy is added to the system. This means that granular materials can stay in a metastable state indefinitely. 
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A schematic "jamming" phase diagram was proposed in Ref. [1^ in order to unite the concepts of jamming in 
many different systems. According to the diagram, jamming in supercooled Uquids occurs at low temperature T and 
high pressure P. Such unified picture suggests several questions, such as whether the dynamics of different systems 
approaching the jammed state are similar. The jamming phase diagram also suggests that the glass and jammed states 
may be related. However, the details of this relationship are not yet fully understood. For instance, a commonly 
measured quantity in granular materials is the probability distribution of normal forces P{F). P{F) develops a peak 
near the jamming transition in granular materials that is well established in experiments and computer simulations. 
It has been suggested that a similar peak in P{F) measured in supercooled liquids signifies the onset of solidity, and 
consequently is a signature of the glass transition. The obvious supposition is: if the link between jamming and the 
glass transition exists, then the mechanism for slowing down of dynamics in supercooled liquids may be related to 
macroscopic structural arrest in granular materials. A peak in PiF) at and below Tg would provide a link in this 
regard. Such a peak was reported in 2-d simulations of several glass- forming liquids [ll . We test the generality of this 
result in the present work. 

As shown in many experiments and computer simulations of granular materials [Ij, [la, [ifl ItTI Il8l | , large force 
inhomogeneities are responsible for structural arrest under shear. If the shear rate is sufficiently low, the granular 
material will be structurally "arrested" or jammed. Jamming in granular materials occurs because the particles form 
"force chains" along the direction in which the stress is applied [13 . This leads us to pose the following question: Do 
force chains exist in supercooled liquids, and if so, are they responsible for the slowing down of dynamics, and related 
to other prominent features of supercooled liquids? 

A salient feature of glass-forming liquids whose origin is still not understood is the spatially heterogeneous nature 
of their dynamics. In particular, the non-exponential character of the relaxation of density correlation functions and 
decoupling of the transport coefficients can be rationalized with the existence of spatially heterogeneous dynamics 
(SHD) or "dynamical heterogeneity" . We refer to a system as dynamically heterogeneous if it is possible to select a 
dynamically distinguishable subset of particles by experiment or com put er simulatio n 1201 . Simulations |2ll 123. 123. 
il m m 113 and experiments [HiimEllilinilHlElS ElM El ES EE Efll have demonstrated the 
cooperative and spatially heterogeneous nature of the liquid dynamics (for reviews of the experimental evidence for 
spatially heterogeneous d yna mics, see, e. g., Refs. [44 \4R l46j ) . 

In Refs. I2, I23, |43j l48L l49l ISOl Isil l53 |. several approaches - including calculation of a displacement-displacement 
correlation function, identification of clusters of mobile particles, and calculation of a four-point time-dependent 
density correlation function - predicted and demonstrated the importance of spatially heterogeneous dynamics in 
supercooled liquids using a rigorous statistical mechanical analysis. In particular, using the four-point time-dependent 
density correlation function formalism, we found dynamical correlation lengths of regions of localized and delocalized 
particles which suggests a picture of fluctuating domains of temporarily localized and delocalized particles, as suggested 
by Stillingcr and Hodgcdon 53] . It has recently been demonstrated in two simulated liquids 54] that these fluctuating 
domains are dynamically facilitated, as predicted by Garrahan and Chandler j55j . Many specific predictions made 
possible through such analyses have now been confirmed in experiments on colloidal suspensions [S^, I13j uM ■ The 
tools from these analyses are thus available to investigate the connection between jamming in granular materials and 
SHD in supercooled liquids. In this work, we will use the four-point time-dependent density correlation function 
formalism to test this connection. 

Here we analyze the forces experienced by particles in a glass-forming liquid in analogy to those in granular materials, 
investigate P{F), and attempt to ascertain a relationship between the forces and spatially heterogeneous dynamics. 
The paper is organized as follows. In Sectional we describe the model and method used to produce our results. In 
Section HTll we give a brief description of SHD and the theoretical framework used to describe SHD. In Section Hvl 
we measure the instantaneous forces between all particle pairs in our glass-forming liquid, and calculate the average 
force and corresponding standard deviation for every T simulated. In Section we divide the set of instantaneous 
pair forces for each configuration into subsets of highest, average and lowest forces at each T. Using our definition 
of localized and replaced particles (defined in Section [fflj, we find the fractions of these particles in each subset of 
instantaneous pair forces. In Section IVIl we define "paths" of highest, average and lowest forces, which we call "force 
chains" , and investigate the temperature dependence of their average mass, average number, and distribution. In 
Section rVlII we investigate the fractions of localized and replaced particles in these "force chains". We conclude with 
the discussion of our results in Section IVIIII 

II. MODEL AND METHOD 

The simulation method we use to generate data for our analyses is molecular dynamics (MD). This is a widely 
used method in the investigation of supercooled liquids and glasses that provides static and dynamic properties for a 
collection of particles described by classical force fields. To perform our simulations we use LAMMPS ^], a parallel 



(T) 


Tq 


0.588 ±0.001 


3500 ± 100 


0.598 ± 0.002 


1900 ± 100 


0.615 ±0.001 


880 ± 50 


0.637 ± 0.001 


370 ± 50 


0.660 ± 0.001 


240 ± 30 


0.689 ± 0.001 


150 ± 30 


0.944 ± 0.001 


16 ±5 


2.004 ± 0.001 


4±1 



TABLE I: Average temperature (T) and relaxation time Ta 



MD code developed by Plimpton. We study a 50/50 binary mixture of particle types "A" and "_B" that interact via 
the Lennard- Jones potential 



Vk/sW = 46^/3 
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This system has been extensively studied by Wahnstrom [53, Schr0der ^3|j a-nd our group |60l l61| . Following 
previous work, we use length parameters aAA = 1, ctbb = 5/6, and ctab = (ctaa + o'bb)/2, and energy parameters 
^AA = ^BB = ^AB = 1- The masses of the particles are chosen to be niA — 2 and ms = 1. We shift the potential and 
truncate it so it vanishes at r = 2.b<TAB- 

We simulate a system of A^ = 8000 particles using periodic boundary conditions in a cubic box of length L — 18.334 
in units of (Taa, which yields a density of p = N/L^ = 1.296 for all state points. We report time in units of 
T = {■mBO'\j^/4:8eAA)^ , length in units of aAA, and temperature, T, in units of eAA/ks, where ks is Boltzmann's 
constant. We simulate state points in the NVE ensemble at temperatures ranging from T — 2.0 to T = 0.001, 
following a path similar to that followed in Refs. |59l kKI l6ll l6^. Additionally, we simulate a series of state points 
in the NPT ensemble for temperatures ranging from 0.1 to 0.6 at zero average pressure. Simulation details can 
be found in Refs. |53. l63|. We estimate the mode coupling temperature Tmct = 0.57 ± 0.01, (the glass transition 
temperature Tg is typically in the range O.QTmct < Tg < O.QTmctISJ]) and the Kauzmann temperature Tq, which 
can be considered a lower bound for the glass transition temperature Tg, is Tq = 0.48 ± 0.02. These quantities are 
found by calculating the structural relaxation time Ta by fitting the secondary relaxation of the coherent intermediate 
scattering function F(qo,t), evaluated at wavevector go corresponding to the maximum peak in the static structure 
factor, to a stretched exponential function F{t) = Acxp{—{t/Ta)^). Table Q] summarizes the T-dependence of Tq for 
the simulations in the NVE ensemble. 



III. BASIC QUANTITIES USED TO MEASURE SPATIALLY HETEROGENEOUS DYNAMICS 



The quantities relevant to SHD that we use here were calculated in Ref. [^ using a theoretical framework based on 
a four-point, time-dependent, density correlation function (74 (r,t). Here, we give definitions of these quantities that 
will be used in later sections of the paper. The quantity g4{r,t) is related to an order parameter Q(t) corresponding 
to the number of "overlapping" particles in a time window t, where the term "overlap" is used to denote a particle 
which was either localized or replaced in a time t. Mathematically, Q{t) is defined as 
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where w(|ri — r2[) is unity for |ri — r2| < a and zero otherwise. We take a — 0.3. The reason for introducing 
an "overlap" function w is to eliminate the vibrational motion of the particles, which is known to be only weakly 
correlated at best (for more details see e.g. Ref. |53|). In Refs. 52, 60], we showed that Q{t) can be decomposed into 
self and distinct components. 
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The self part, Qs{t), measures the number of particles that move less than a distance a in a time interval t; we call 
these "localized" particles. The distinct part, Qnit) measures the number of particles replaced within a radius a by 
another particle in time t; we call these "replaced" particles. 

We also consider "delocalized" particles, that is, particles that in a time t moved more than a distance a from their 
original location. As was pointed out in Ref. [63], substituting 1 — w for it; in Eq. ((2Jl gives the delocalized order 
parameter Qohii) — N ~ Qs{t)- 

The fluctuations, XiW^ in Q{t) niay be defined as 

x4{t)^^[{Qm-{Qit)n (4) 

Following the scheme of decomposing Q{t), x^it) can be decomposed into self Xssit), distinct Xooit), and cross 
Xsoit) terms. xss{t) — XDL{i) is the susceptibility arising from fluctuations in the number of locahzed particles, 
XDoit) is the susceptibility arising from fluctuations in the number of particles that are replaced by a neighboring 
particle, xsoit) represents cross fluctuations between the number of localized and replaced particles, and xdl is the 
susceptibility arising from fluctuations in the number of delocalized particles. 

Xiit) (and its terms) measures the correlated motion between pairs of particles, calculated equivalently from fluctu- 
ations in the number of overlaps or from the four-point correlation function itself. As shown in Ref. 52] , the behavior 
of X4(i) demonstrates that correlations are time dependent, with a maximum at a time i™'*''. The characteristic times 
t^^^ and Ta (defined in the previous section) have similar T-dependence, indicating that the correlations measured 
by x^it) 3'i'6 most pronounced in the structural a-relaxation regime. 

IV. AVERAGE FORCE AND INSTANTANEOUS FORCE DISTRIBUTION FUNCTION P{F) 

Because we are studying a dense liquid described by a relatively short range pair potential, at an instant in time 
each particle experiences a force due to the presence of neighboring particles. Following Ref. |3], we calculate the 
average force (F) between every pair of neighboring particles at a given instant of time. In Figure ^ we show (F) 
vs T calculated from at least iVconf = 1000 independent NVE configurations at each T. We see that both {F) and 

the standard deviation cr^^) — ( jj t^^^ —pr ^^ ( {F)c — {P})^ ) , where {F)c is the average force for a particular 

configuration (c), decrease with decreasing T. [We also find that the standard deviation <y(F) calculated within each 
configuration has the same temperature dependence (not shown).] A decrease is expected since the average pressure 
(see Table 1 in Ref [sj) in our constant volume simulation decreases as T decreases. Large a(^p^ at higher T is also 
expected because of the larger fluctuations in pressure at higher temperatures. (F) is calculated in the following 
manner. First, the average force of a particular configuration c is found from {F)c — (X^ij Pij)/^ijj where Nij is the 
number of pairs i and j for which the force is nonzero, and particles i and j belong to that particular configuration. An 
ensemble average force (F) is calculated as an equally weighted average over all {F)c. This can be an important point 
when there are large force fluctuations that result in substantially different numbers of positive forces for different 
configurations. We have, however, checked an alternative method to calculate average force by calculating the average 
force from all particle pairs making no distinction between configurations, and obtained the same answer, indicating 
that the system is self-averaging. In Ref. |65|] the authors showed that the presence of non-self-averaging alters the 
force distribution function P{F), which we now consider. 

In foams and granular materials, P{F) is measured as a distribution of interparticle normal forces |6a,|63. In our 
system, there is no strict definition of point contact between Lennard-Jones(LJ) particles, and therefore the normal 
force between particles in not well-defined. We define two LJ particles to be in contact if they interact with a positive 
(repulsive) force, as it was done in Ref. 1]. Commonly, granular materials are modeled as systems of particles with 
only repulsive interactions when grains are in contact. Therefore, we investigate P{F) only for the subset of positive 
(repulsive) F. Figure [21 shows the T-dependence of the force distribution function P{F), calculated as a normalized 
histogram of all instantaneous positive forces between particle pairs. The behavior of P{F) is similar to the one 
observed in Ref.jj for their LJ system above Tg. We observe an exponential tail at high forces and an increase in 
the curvature of P{F) at small forces as T decreases. An explanation of the origin of the long exponential tail in 
P{F) was proposed in Ref. [l| and is related to the fact that large F behavior can be obtained from the small r 
behavior of g{r) « y(r)exp[— V(r)/kbT], and P{F)dF oc q (r)dr. The exponential tail in the force distribution of 
granular materials was suggested in the g-model in Ref. |68j ] to be a consequence of a force randomization throughout 
the packing. This randomization has an effect analogous to the collisions in an ideal gas ,68. .69.] . 
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FIG. 1: T-dependence of the average force (F). The inset shows the T-dependence of the standard deviation cr^^^. 
bars on (F) and a/^p^ are smaller than the symbol size. 
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To further compare P{F) with that in granular materials, we calculate the force distribution function scaled with 
the average force {F)c for a particular configuration, P{F/{F)c)- The T dependence of P{F/{F)c) vs F/{F)c is 
shown in Figure O Wc include P{F/{F)c) for several temperatures below Tq = 0.48 for comparison with the force 
distribution found in Rcf. 1] below Tg. As in Figure[21 the probability of observing forces substantially larger than {F) 
decays exponentially, and the curvature around (F) decreases as T decreases (See FigureOlcaption). The resemblance 
of Figure 121 and Figure to similar figures in Ref. [l| again shows there is no significant distinction between averaging 
pair forces within a configuration and globally (across configurations), confirming that large force fiuctuations, like 
those found in Ref. |65| for an out-of-equilibrium system, are not present here. 

We note that the reason for using two methods for the calculation of P{F) is the fact that glass-forming liquids, such 
as our model, are not in equilibrium on all time scales. Long lived clusters of immobile particles have been reported 
to persist for times that are orders of magnitude longer that the structural relaxation time r^ |70|| . It is intuitive 
that if those clusters exist for a long time, forces between mobile particles could have large fluctuations in order to 
initiate structural relaxation. Therefore, it is possible that large force fluctuations will be present in supercooled 
liquids even at structural relaxation times. If this was the case, then one could observe non self-averaging of P{F). 
We demonstrate that this is not the case by comparing two different methods of calculating P{F). 

In the experiments and simulations of granular materials and foams, it has been suggested that a signature of 
jamming is the development of a small peak or plateau around the average force. O'Hern et al. suggested that P{F) 
will develop a peak if the first peak of g{r) is "sufficiently high and narrow" [71| which they associate with solid-like 
behavior of the system and a possible signature of the glass transition. Indeed we observe the sharpening of the first 
peak in g{r) (not shown here, see e. g. Ref. |73|) in our system as T decreases, but the peak in P{F) is absent for 
T < Tq. Recall that the jamming transition in granular materials is equivalent to the macroscopic structural arrest of 
the system. This would correspond to a glass transition in supercooled liquids. Thus for T > Tg, we would not expect 
to observe the peak or plateau in P{F/{F)c) (see Figurel^J. In order to investigate the behavior of P{F/ {F)c) in the 
glass, we calculate P{F/{F)c) at several T below Tq. These state points are obtained by quenching a configuration 
previously equilibrated at high temperature, in our case T = 10.0, to a desired T, at a non-zero pressure. We do 
observe that the curvature of P{F/{F)c) around (F) decreases as T decreases, and that the slope of logP{F/{F)c) 
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FIG. 2: T dependence of P{F) vs F. The inset shows behavior of P{F) for smaU positive forces. 



for 10 < (F) < 60 changes significantly for T = 0.001. We observe that the peak or plateau is absent in our system 
approaching Tmct from above, and that the plateau develops at much lower temperatures (T = 0.1) compared to 
To. This demonstrates that development of a peak or plateau is not a necessary condition for a glass transition or a 
solid- like behavior in dense glass- forming systems. 

In the simulations of Ref. [l|, a barostat was used to maintain an average zero pressure. To test if the peak 
development in P{F / {F)c) depends on the choice of thermodynamic ensemble, we performed simulations of our 
3-D binary LJ system at several temperatures in the NPT ensemble, and calculated the corresponding P{F/{F)c). 
Figure 01 shows the T-dependence of P{F/{F)c) for several NPT state points. The pressure is kept at zero for all 
temperatures. We again observe a change in the curvature in P{F), but only at very low temperatures below Tg, 
consistent with what our findings for the NVE system. This again shows that a peak in P{F/{F)c) is not necessary 
for solid-like behavior of our system. Similar findings have been reported by Reichman and Sastry for their model 
supercooled liquid J7^. 

It is interesting to note that both the NVE and NPT force distributions P{F/{F)c) have a value of F/{F)c where 
the distributions cross for almost all T, indicating an isosbestic point. This value is at approximately F/{F)c — 10 and 
55 for the NVE and NPT systems, respectively. The reason for such behavior is not obvious, and will be investigated 
in future work. Also note that we tested the self-averaging of P{F/{F)c) for all state points for which the system is 
a glass. We find no evidence for non self averaging of P{F/{F)c) at those points. 

Despite the fact that the peak or plateau in P{F) is observed in many granular materials and some supercooled 
liquids below Tg, an open question is its connection to the development of the yield stress in those materials. Several 
speculations have been made in Refs. [l|, l7lj | in order to explain the jamming scenario and its possible connection 
to the glass transition. They theorize that "systems jam when there are enough particles in a force chain network 
to support stress over the time scale of the measurement" which would imply that force chains observed in granular 
packing may also be important to the glass transition. Further, since force chains do not couple strongly to density 
fluctuations they may be linked to local dynamical heterogeneities near Tg. To test this speculation we seek to find 
a link between SHD and particles that belong to the different regions of the force distribution function, namely the 
particle pairs that interact with forces that belong to the exponential tail and average force in the force distribution 
function. 
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FIG. 3: T-dependence of the force distribution P(F/(F)c) vs F/{F)c- Forces that are much larger than (F) decay exponentially. 
The curvature of P{F/{F)c) around F/{F)c = 1 does not change sign but its value decreases as T decreases. The inset shows 
a close-up of P{F)/{F)c for forces near {F)c. 



SUBSETS OF INSTANTANEOUS FORCES 



In Ref 5z|, we established a criterion to find localized and replaced particles which we use here to test the relationship 
between particle mobilities and instantaneous pair forces. A brief description of how to identify localized, replaced 
and delocalized particles can be found in Section UTTI For every two configurations (one at t = and the other at t) we 
find the number of localized {Qs(t)) and replaced {Qoii)) particles and calculate all instantaneous pair forces in the 
configuration at time t. We sort all instantaneous pair forces according to their values and find subsets of particles 
that interact with pair forces that fall within certain percentages of the lowest, average and highest pair forces. The 
considered percentages of these forces are indicated in Figure 

We define a fraction, $ioc ($rep), as the ratio of the number of localized (replaced) particles that are associated 
with the given subset of pair forces to the total number of localized (replaced) particles, Qs (Qd)- 'i'loc and ^j-cp thus 
relate mobility and instantaneous forces. Figure |S1 shows the time dependence of $ioc and ^icp in the given subsets 
of pair forces at T = 0.60. 

For the limiting case where all localized and replaced particles belong to the given subset of pair forces, <I>ioc = 1 
and <I>rop = 1. This is true for subsets larger than « 5% of the pair forces. These fractions are large enough to include 
nearly all of the system's particles since there are many more pair forces (w 225000) than particles (8000). For lower 
percentages of the highest pair forces, ^loc decreases in time from the value expected for the bulk (at t=0 when all 
particles are localized and <f>ioc is simply the probability of finding any particle in the given percentage of pair forces) . 
This suggests that because particles that have been localized for long times are less likely to experience high relative 
forces, if the reason for their long localization is due to being at a local energy minimum. $iop decreases in time 
toward the bulk value (at later times), but this decrease starts on much shorter time scales than $ioc, possibly because 
particles that have been replaced (mobile particles) experience higher relative forces at times when they escape from 
their cages. In the case of low pair forces, the results for $ioc and <i>rep should be taken with caution because there 
are two ranges of distance at which particles can experience the lowest force (the tail and well of the potential). The 
ambiguity in low forces is seen in the inconsistent behavior of ^loc and 'I'rcp- This may mask any clear interpretation 
of the meaning of $ioc and 'I'rop for the lowest forces. In the case of average pair forces, the results are reversed from 
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FIG. 4: T-dependence of the force distribution P(F/(F)c) vs F/{F)c at several state points in the NPT ensemble at P = 0. 
The inset shows a close-up oi P(F)/{F)c for forces near {F)c- 

those of the high pair forces. This means that a larger fraction of localized particles tend experience average forces 
at later times, and replaced particles tend to approach the bulk value of the average force. 

This connection between mobility and instantaneous forces leads us to investigate the spatial correlations of pairwise 
forces in supercooled liquids. 



VI. GENERALIZED FORCE CHAINS IN SUPERCOOLED LIQUIDS 



We now investigate the spatial correlations among the instantaneous pairwise forces. Spatially correlated forces 
can be viewed as networked "paths" of highest, average or lowest forces in the system. Highest-force paths would 
correspond to the force chains in granular materials, and we refer to them as force chains. 

According to Ref . [TJI , a force chain in a granular material is a "linear string of rigid particles in point contact" . 
A force chain in a granular material can support a load along its axis. Since we do not apply a load or stress to our 
systeni, if such structures exist, they would be induced by lowering the temperature. Inspired by the work of Makse 
et al. upt . we look for these chains by searching the "paths" that may be found following instantaneous forces. In 
Ref. [73 , the authors found force chains by starting from a sphere at the top of the box of simulated spherical grains 
and following the path of maximum contact force at every grain. They observed a force-bearing network that was 
concentrated in a few percolating chains. 

In our case, we define force chains such that, for each particle, we find the two neighboring particles that exert on 
that particle the highest force. This set of three particles constitute what we call a "trivial chain" , since it is always 
formed by construction. These trivial chains are schematically shown in Fig. |3 Later, chains that are defined using 
low and average forces are also calculated in a similar manner as chains defined using the highest forces. Longer 
chains are formed from trivial chains that share two members (Figure [TJ. Force chains from the highest forces would 
correspond to force chains of point contacts in granular materials. We investigate low and average force chains for the 
sake of completeness and the fact that we do not have a theory that would a priori suggest whether jamming occurs 
because of high forces or, e. g., average forces. 

Figure IHl shows a chain of highest forces containing 26 particles in a single configuration at T — 0.60. Figure 
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FIG. 5: Time dependence of "I>ioc (upper three panels) and $rcp (lower three panels) in the subsets of low, average and high 
pair forces at T = 0.60. $rcp is shown for times when Qd becomes nonzero. Labels on the y axes for average and high forces 
are omitted because of the clarity. 





FIG. 6: Schematic picture of trivial force chains. Each particle is connected to two neighbors with which it interacts with the 
highest forces. Fji and Fjk are the largest forces on particle j. A similar picture can be drawn for interactions with low or 
"closest to average" forces. 
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FIG. 7: Schematic picture of formation of non-trivial chains. This is an example of a chain of mass six. In this example four 
trivial chains combine to form a chain of mass six. 



shows all highest force chains at T = 0.60 at one snapshot. We refer to chain mass n as the number of particles in 
that chain. A'chains represents the number of non-trivial chains for a particular configuration. 

The first trend we investigate is the T-dependence of the average mass (n) of the force chains and the average 
number of non-trivial chains (A'chains}- These quantities are shown in Figure^] (n) and (A'chains} are averages over 
the same configurations we used to calculate the force distribution functions, P{F), in Section Hvl It is evident from 
Figure ITUI that (n) and (iVchains) have at most a very weak temperature dependence. We observe that chains with 
high forces are on average slightly longer than chains with average forces. We also observe that the average number 
of chains calculated for the average force subsets is slightly larger than the average number of chains calculated for 
low and high forces. 

In order to look for any spanning chains (similar to Ref. |75|), we calculate the radius of gyration and end-to-end 
distance of the force chains, common quantities used to describe the size and shape of polymers [76J. We find that 
these quantities are also largely independent of temperature (not shown). Therefore, the geometry of force chains 
does not change significantly with T. We also do not find any spanning chains that would correspond to the case 
where the end-to-end distance is equal to the size of the box or where the radius of gyration equals half of the box 
size. 

In order to better understand the difference between chains that are defined using different subsets of forces, we 
investigate the distribution function of chain mass P{n) for chains defined using low, average and high forces. P{n) 
is the probability of finding chains of mass n using the given subset of forces. P{n) is calculated as a normalized 
histogram of chain masses for each configuration and averaged over all configurations. Figure [TTl shows P{n) vs n for 
chains defined for each subset. We see from P{n) that chains of high forces are, on average, longer than chains of 



11 




FIG. 8: 26-particle chain ai T — 0.60. The particles in the chain are shown with radius 0.5cr. This is the longest chain for this 
particular snapshot. This chain does not span the box. Dots represent the rest of the 8000 particles in the system. 
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FIG. 9: Snapshot of all force chains at T=0.60. Different colors correspond to force chains with different masses. 



average and low forces, as seen in Figure [TUT a). Shorter chains defined using average forces may be the consequence of 
more branching of these chains compared to the chains composed of the highest and lowest forces, which is consistent 
with the fact that the number of force chains for the average forces is the largest. 

To give a theoretical prediction of P{n) for comparison, consider the case where the forces between particles are 
randomly assigned. Assume that there are k nearest neighbors for a single particle X. If a particle, X, shares one of 
its highest pair forces with a neighboring particle, A, the probability that the interaction between A and X is one of 

A's two highest pairwise forces can be approximated as ( f ) . Assuming that the chain continues to grow away from 

its origin and the likelihood of it looping back on itself is small, the probability of adding ni additional members on 
one end of a trivial chain {n = 3) f Figure [T^ is 



P(ni) 



(-!)■ 



(5) 



Since there are two ends to the chain and the trivial chain always has three members, the probability for a total chain 
length 71 = ni + 712 + 3 is 



^(")=("-2)©"" 
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FIG. 10: T dependence of the (a) average mass n of the force chains and (b) average number (A'^chains) of non-trivial force 
chains for low, average and high pair forces. The error bars are smaller than the symbol size. 



where the linear factor n — 2 accounts for all combinations of rii and n2 that sum to n — 3. Figure lTTI shows predictions 
of P{n) for various values of the nearest neighbor number parameter, k. From the integration of g(r), the actual 
number of neighbors in the first neighbor shell is approximately 13. As seen from the figure, the measured P{n) 
decays much more slowly than the prediction for entirely randomly assigned forces. Instead, the observations are well 
bounded by k values of three and four. This gives some measure of the variability in the local force distribution. In 
other words, the more the local environment of one particle varies from that of its neighbor, the larger the value of 
k that is found. We also note the possible connection between the force chains found here and the idea of rigidity 
percolation by Phillips [77| and Thorpe [78j. In rigidity percolation, the mean coordination number m = 2.4 represents 
a threshold below which a network glass is easily deformed. The fact that our nearest neighbor parameter k is close 
to three suggests the possibility that the network of bonds in network glasses and the network of force chains in the 
LJ mixture studied here share similar properties, and one could use ideas developed in rigidity percolation theory to 
study force chains in supercooled, non network-forming liquids. 

Figurein^shows the T-dependence of P{n) for chains in low, average and high force subsets. We find that P{n) can 
be generally fitted well by a functional form P{n) — aiexp[— a2n], where ai and 02 are fitting parameters. We note 
that observed exponential behavior of P{n) is analogous to that reported for equilibrium polymerization 79] of linear 
polymer chains, in which the bonds between monomers break and recombine at random points along the backbone 
of the chain. This picture is exactly what happens to the force chain network in our system. Furthermore, the size 
distribution of strings of cooperative rearranging particles is exponential in all studies of strings performed thus far. 
In the case of chains defined using average forces we note a slight T-dependence of P{n). It appears that these chains 
are shorter at lower temperatures, but the distribution is still exponential. In the case of the chains defined using the 
smallest forces, chains with mass n = 4 and n = 5 do not fit well with the exponential function (note the slight bend 
in the curve in Figure ITST a)). 

We note that the force chains, as constructed, may be highly susceptible to thermal fluctuations. More subtle trends 
may be detectable by using methods to minimize these fluctuations such as quenching to the inherent structure or 
time averaging over time scales similar to Tq, so as to filter out frequencies higher than those that act over the periods 
of structural rearrangement. 

We conclude this section by summarizing that we do not find a strong T dependence of the chains defined above. 
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FIG. 11: Comparison of the measured mass distributions of force chains for low, average and high force at T = 0.60 and 
theoretical predictions for different nearest neighbor parameter, k, values. 

This is surprising because if the highest forces (or any force chains) are related to slowing down of dynamics in 
supercooled liquid, one might expect to see a temperature dependence pattern in the chain properties. 

VII. FRACTION OF LOCALIZED AND REPLACED PARTICLES IN FORCE CHAINS 



To look for a connection between force chains and mobility, we calculate the fraction of localized and replaced 
particles in each chain for all subsets of forces for which force chains are defined. This fraction is calculated in the 
following manner. For each configuration at t, we identify the force chains and localized and replaced particles with 
respect to an initial configuration at t — 0, and for each chain of mass n in the configuration at t we count the number 
of localized and replaced particles and divide by n. After this, we average those fractions over all equally spaced 
configurations. Figure [HI shows the fraction of localized and replaced particles in the chains of high, average and low 
force at tf'''' (defined for each T in Ref ^52], and Section Uni) at T = 0.60, T = 0.66 and T = 0.94. We examine the 
configurations at the particular time t^'^^ because, as explained in Section ITm SHD as measured by X4(i) and ^4{t) is 
most pronounced then. If SHD in supercooled liquids and granular materials share common mechanisms, one might 
expect the effects of jamming to be most detectable at this time. The fraction of localized and replaced particles is 
shown in Figure^] for each T. The first column in Figure [T^ represents <&chain at T = 0.60 for low, average and high 
force chains. The second and third column contain data for <i>chain at T = 0.66 and T = 0.94 respectively. 

From Figure ^^ we see that $chain is essentially constant for n < 20 for both localized and replaced particles. 
For chain masses greater then 20, ^chain becomes noisy because of poor statistics. For chains with n < 20, $chain 
has values slightly higher than the fraction of localized particles in the bulk, Qs {Qs is indicated by solid lines in 
Figure EJ. This means that localized particles are more likely to be found in chains with intermediate mass than 
would be expected from the fraction of localized particles, (Qs), in the system. $chain has values equal to or slightly 
lower than the fraction of replaced particles in the bulk {Qo is also indicated by solid lines in Figure^J. This means 
that replaced particles are equally or less likely to be found in the chains of intermediate mass based on their bulk 
population. These behaviors are basically universal regardless of the subset of forces in which the force chains are 
defined. Therefore, we cannot make a simple connection with mobility, and force chains as defined here probably do 
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FIG. 12: Example of chain length prediction. 

not play the same role as those found in granular materials. 

The population of localized and replaced particles in the force chains can be explained if we suppose that localized 
particles reside in a local environment that is less susceptible to force perturbations than delocalized particles. Imagine 
an " intrinsic" high force chain that would be associated with the inherent structure and introduce a perturbation near 
the chain that propagates perpendicular to it. This perturbation can temporarily change which neighbors of a given 
particle in the chain interact with it via the highest relative force, thereby breaking the chain. Such perturbations are 
the manifestation of temperature in the system. Since the localized particles are more likely to be found in non-trivial 
chains, one might speculate that something about their local environment makes them resistant to these perturbations 
at least with respect to the relative forces exerted by their neighbors. 

VIII. DISCUSSION 

In this paper we examined the relationship between jamming in granular materials and SHD in supercooled liquids. 
We calculated the instantaneous force distribution function P{F) in a model glass-forming liquid, and we did not find 
a peak in this quantity at any T whether above or below Tg, in contrast to the model supercooled liquids studied in 
Ref. [l|. We also found a possible connection between instantaneous force magnitude and long term mobility. We 
defined force chains in our model glass-forming liquid, and based on our results from Section IVll we found force chains 
much longer than those that might be expected for randomly assigned forces. These force chains probably do not 
play the same role in this supercooled liquid as in granular materials because they do not show a strong temperature 
dependence. 

Force chains may, however, indicate a difference in the evolution of the local environment of particles with different 
mobilities which may be connected with cooperative and string like motion found in model supercooled liquids close to 
Tg. Microscopic details of local particles dynamics and the mechanism by which particles move along string- like paths 
is studied in Ref. Q in a glass-forming Dzugotov liquid. The authors show that simultaneous motion of individual 
particles along the string depends on the length of the string, and that for shorter strings the motions is highly 
coherent and for longer strings motion is coherent only within short segments containing as many as seven particles 
(micro strings). We note the similarity between the anatomy of strings and force chains studied here by comparing 
a snapshot of a string in Figure 10 of Ref. Q and Figure |H1 of this paper (remember that the criterion for a string 
is completely different from the criterion for a force chain) . We also note that both the distribution of mass of force 
chains and distribution of string length appear to be exponential. Although the mass of the force chains docs not grow 
as a function of T, in contrast to string length, these distributions are, to our knowledge, the only two distributions 
of quantities associated with force and particle mobility measured in supercooled liquids that obey exponential laws. 
This leads us to pose certain questions: What is the relationship between force chains and strings? Is it possible 
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FIG. 13: T dependence of the chain mass distribution P{n) for the (a) low, (b) average and (c) high force. Force chains defined 
using the average force subset show a shght T-dependence. 

that force chain networks provide a characteristic structural feature along which strings may occur? Could frequent 
defects in the force chain network due to breakage and recombination of trivial force chains be associated with the 
local liquid excitations proposed by Garrahan and Chandler |53 , and provide a path for string like motion? Answers 
to those questions may help us to understand the origin of SHD. Further work on the relationship between strings 
and force chains is ongoing. 
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